// This code replicates tables and figures in the paper titled "Estimating 
// Consumer Substitution between New and Used Passenger Vehicles" published in 
// the Journal of the Association of Environmental and Resource Economists 

// ========================================================================== //
// ================================ Tables ================================== //
// ========================================================================== //
clear all

cd "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021"
use "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021\tables1_4data.dta"


// ================================ Table 1 ================================= //
summarize sales 
summarize trans_price cpm hp_weight footprint new_2 used_2 AWD Sedan SUV Hybrid PEV [w=sales]



// ================================ Table 2 ================================= //

// OLS - No fixed effects
xi: reg lnsales lnprice cpm hp_weight footprint my_avg, vce(cluster model)

// OLS - With fixed effects
xi: reg lnsales lnprice cpm hp_weight footprint my_avg i.drive_type_w i.drive_type_w  i.fuel_type_w i.body_style_w i.cyl_w, vce(cluster model)

// IV - No fixed effects
xi: ivregress 2sls lnsales (lnprice = cpm_diff_in cpm_diff_out hp_weight_diff_in hp_weight_diff_out  footprint_diff_in footprint_diff_out) cpm hp_weight footprint my_avg, vce(cluster model)

// IV - With fixed effects
xi: ivregress 2sls lnsales (lnprice = cpm_diff_in cpm_diff_out hp_weight_diff_in hp_weight_diff_out  footprint_diff_in footprint_diff_out) cpm hp_weight footprint my_avg i.drive_type_w i.drive_type_w  i.fuel_type_w i.body_style_w i.cyl_w, vce(cluster model)

local epsilon = _b[lnprice]


// Market Price Elasticity
egen total_sales = sum(sales)

gen market_elasj = delta*`epsilon'
replace market_elasj = market_elasj*sales/total_sales
egen market_elas = sum(market_elasj)

summarize market_elas



// =============================== Table 4 ================================== //
gen epsilon_hat = `epsilon'

// Own-price elasticity = -3, diversion fraction = 50% of baseline
gen epsilon1 = -3
gen delta1 = 0.5*delta

gen market_elasj1 = delta1*epsilon1
replace market_elasj1 = market_elasj1*sales/total_sales
egen market_elas1 = sum(market_elasj1)

summarize market_elas1


// Own-price elasticity = -3, diversion fraction = baseline
gen epsilon2 = -3
gen delta2 = delta

gen market_elasj2 = delta2*epsilon2
replace market_elasj2 = market_elasj2*sales/total_sales
egen market_elas2 = sum(market_elasj2)

summarize market_elas2


// Own-price elasticity = -3, diversion fraction = 2*baseline
gen epsilon3 = -3
gen delta3 = 2*delta

gen market_elasj3 = delta3*epsilon3
replace market_elasj3 = market_elasj3*sales/total_sales
egen market_elas3 = sum(market_elasj3)

summarize market_elas3



// Own-price elasticity = benchmark, diversion fraction = 50% of baseline
gen epsilon4 = epsilon_hat
gen delta4 = 0.5*delta

gen market_elasj4 = delta4*epsilon4
replace market_elasj4 = market_elasj4*sales/total_sales
egen market_elas4 = sum(market_elasj4)

summarize market_elas4


// Own-price elasticity = benchmark, diversion fraction = baseline
gen epsilon5 = epsilon_hat
gen delta5 = delta

gen market_elasj5 = delta5*epsilon5
replace market_elasj5 = market_elasj5*sales/total_sales
egen market_elas5 = sum(market_elasj5)

summarize market_elas5


// Own-price elasticity = benchmark, diversion fraction = 2*baseline
gen epsilon6 = epsilon_hat
gen delta6 = 2*delta

gen market_elasj6 = delta6*epsilon6
replace market_elasj6 = market_elasj6*sales/total_sales
egen market_elas6 = sum(market_elasj6)

summarize market_elas6



// Own-price elasticity = -6, diversion fraction = 50% of baseline
gen epsilon7 = -6
gen delta7 = 0.5*delta

gen market_elasj7 = delta7*epsilon7
replace market_elasj7 = market_elasj7*sales/total_sales
egen market_elas7 = sum(market_elasj7)

summarize market_elas7


// Own-price elasticity = -6, diversion fraction = baseline
gen epsilon8 = -6
gen delta8 = delta

gen market_elasj8 = delta8*epsilon8
replace market_elasj8 = market_elasj8*sales/total_sales
egen market_elas8 = sum(market_elasj8)

summarize market_elas8


// Own-price elasticity = -6, diversion fraction = 2*baseline
gen epsilon9 = -6
gen delta9 = 2*delta

gen market_elasj9 = delta9*epsilon9
replace market_elasj9 = market_elasj9*sales/total_sales
egen market_elas9 = sum(market_elasj9)

summarize market_elas9



// =============================== Table 5 ================================== //
clear all

// CAFE, 3% discount rate, elasticity = -1 
insheet using "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021\annual_societal_costs_report_cafe_-1.csv"

// 3% discount rate
keep if discrate >= 0.025 & discrate <= 0.035

// Keep baseline and 1.5% scenario
keep if scenario == 0 | scenario == 4

// Keep total costs
keep if regclass == "TOTAL"

// Sum fatality and nonfatal crash costs
gen accidentcosts = fatalitycosts + nonfatalcrashcosts

drop if modelyear > 2029

collapse (sum) accidentcosts netsocialbenefits, by(scenario)

gen diff_accidentcosts = accidentcosts - accidentcosts[_n-1]
gen diff_netsocialbenefits = netsocialbenefits - netsocialbenefits[_n-1]

summarize diff_accidentcosts diff_netsocialbenefits

clear all


// CAFE, 3% discount rate, elasticity = -0.37 
insheet using "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021\annual_societal_costs_report_cafe_-0.37.csv"

// 3% discount rate
keep if discrate >= 0.025 & discrate <= 0.035

// Keep baseline and 1.5% scenario
keep if scenario == 0 | scenario == 4

// Keep total costs
keep if regclass == "TOTAL"

// Sum fatality and nonfatal crash costs
gen accidentcosts = fatalitycosts + nonfatalcrashcosts

drop if modelyear > 2029

collapse (sum) accidentcosts netsocialbenefits, by(scenario)

gen diff_accidentcosts = accidentcosts - accidentcosts[_n-1]
gen diff_netsocialbenefits = netsocialbenefits - netsocialbenefits[_n-1]

summarize diff_accidentcosts diff_netsocialbenefits

clear all


// CAFE, 7% discount rate, elasticity = -1 
insheet using "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021\annual_societal_costs_report_cafe_-1.csv"

// 7% discount rate
keep if discrate >= 0.065 & discrate <= 0.075

// Keep baseline and 1.5% scenario
keep if scenario == 0 | scenario == 4

// Keep total costs
keep if regclass == "TOTAL"

// Sum fatality and nonfatal crash costs
gen accidentcosts = fatalitycosts + nonfatalcrashcosts

drop if modelyear > 2029

collapse (sum) accidentcosts netsocialbenefits, by(scenario)

gen diff_accidentcosts = accidentcosts - accidentcosts[_n-1]
gen diff_netsocialbenefits = netsocialbenefits - netsocialbenefits[_n-1]

summarize diff_accidentcosts diff_netsocialbenefits

clear all


// CAFE, 7% discount rate, elasticity = -0.37 
insheet using "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021\annual_societal_costs_report_cafe_-0.37.csv"

// 7% discount rate
keep if discrate >= 0.065 & discrate <= 0.075

// Keep baseline and 1.5% scenario
keep if scenario == 0 | scenario == 4

// Keep total costs
keep if regclass == "TOTAL"

// Sum fatality and nonfatal crash costs
gen accidentcosts = fatalitycosts + nonfatalcrashcosts

drop if modelyear > 2029

collapse (sum) accidentcosts netsocialbenefits, by(scenario)

gen diff_accidentcosts = accidentcosts - accidentcosts[_n-1]
gen diff_netsocialbenefits = netsocialbenefits - netsocialbenefits[_n-1]

summarize diff_accidentcosts diff_netsocialbenefits

clear all


// EPA, 3% discount rate, elasticity = -1 
insheet using "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021\annual_societal_costs_report_epa_-1.csv"

// 3% discount rate
keep if discrate >= 0.025 & discrate <= 0.035

// Keep baseline and 1.5% scenario
keep if scenario == 0 | scenario == 4

// Keep total costs
keep if regclass == "TOTAL"

// Sum fatality and nonfatal crash costs
gen accidentcosts = fatalitycosts + nonfatalcrashcosts

drop if modelyear > 2029

collapse (sum) accidentcosts netsocialbenefits, by(scenario)

gen diff_accidentcosts = accidentcosts - accidentcosts[_n-1]
gen diff_netsocialbenefits = netsocialbenefits - netsocialbenefits[_n-1]

summarize diff_accidentcosts diff_netsocialbenefits

clear all


// EPA, 3% discount rate, elasticity = -0.37 
insheet using "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021\annual_societal_costs_report_epa_-0.37.csv"

// 3% discount rate
keep if discrate >= 0.025 & discrate <= 0.035

// Keep baseline and 1.5% scenario
keep if scenario == 0 | scenario == 4

// Keep total costs
keep if regclass == "TOTAL"

// Sum fatality and nonfatal crash costs
gen accidentcosts = fatalitycosts + nonfatalcrashcosts

drop if modelyear > 2029

collapse (sum) accidentcosts netsocialbenefits, by(scenario)

gen diff_accidentcosts = accidentcosts - accidentcosts[_n-1]
gen diff_netsocialbenefits = netsocialbenefits - netsocialbenefits[_n-1]

summarize diff_accidentcosts diff_netsocialbenefits

clear all


// EPA, 7% discount rate, elasticity = -1 
insheet using "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021\annual_societal_costs_report_epa_-1.csv"

// 7% discount rate
keep if discrate >= 0.065 & discrate <= 0.075

// Keep baseline and 1.5% scenario
keep if scenario == 0 | scenario == 4

// Keep total costs
keep if regclass == "TOTAL"

// Sum fatality and nonfatal crash costs
gen accidentcosts = fatalitycosts + nonfatalcrashcosts

drop if modelyear > 2029

collapse (sum) accidentcosts netsocialbenefits, by(scenario)

gen diff_accidentcosts = accidentcosts - accidentcosts[_n-1]
gen diff_netsocialbenefits = netsocialbenefits - netsocialbenefits[_n-1]

summarize diff_accidentcosts diff_netsocialbenefits

clear all


// EPA, 7% discount rate, elasticity = -0.37 
insheet using "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021\annual_societal_costs_report_epa_-0.37.csv"

// 7% discount rate
keep if discrate >= 0.065 & discrate <= 0.075

// Keep baseline and 1.5% scenario
keep if scenario == 0 | scenario == 4

// Keep total costs
keep if regclass == "TOTAL"

// Sum fatality and nonfatal crash costs
gen accidentcosts = fatalitycosts + nonfatalcrashcosts

drop if modelyear > 2029

collapse (sum) accidentcosts netsocialbenefits, by(scenario)

gen diff_accidentcosts = accidentcosts - accidentcosts[_n-1]
gen diff_netsocialbenefits = netsocialbenefits - netsocialbenefits[_n-1]

summarize diff_accidentcosts diff_netsocialbenefits



// ==================== Appendix Table A.2 Column 2 ========================= //
clear all

cd "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021"
use "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021\tablea2column2data.dta"


// IV - With fixed effects
xi: ivregress 2sls lnsales (lnprice = cpm_diff_in cpm_diff_out hp_weight_diff_in hp_weight_diff_out  footprint_diff_in footprint_diff_out) cpm hp_weight footprint my_avg flexible hybrid phev awd SUV pickup cyl3 cyl4 cyl5 cyl6 cyl8 cyl10, vce(cluster model)

local epsilon = _b[lnprice]


// Market Price Elasticity
egen total_sales = sum(sales)

gen market_elasj = delta*`epsilon'
replace market_elasj = market_elasj*sales/total_sales
egen market_elas = sum(market_elasj)

summarize market_elas



// ===================== Appendix Table A.2 Column 3 ======================== //
clear all

cd "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021"
use "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021\tablea2column3data.dta"

// IV - With fixed effects
xi: ivregress 2sls lnsales (lnprice = cpm_diff_in cpm_diff_out hp_weight_diff_in hp_weight_diff_out  footprint_diff_in footprint_diff_out) cpm hp_weight footprint i.fuel_type_w i.drive_type_w i.body_style_w i.cyl_w, vce(cluster model)

local epsilon = _b[lnprice]


//  Market Price Elasticity 
egen total_sales = sum(sales)

gen market_elasj = delta*`epsilon'
replace market_elasj = market_elasj*sales/total_sales
egen market_elas = sum(market_elasj)

summarize market_elas





// ========================================================================== //
// =============================== Figures ================================== //
// ========================================================================== //
clear all


// ============================== Figure 1 ================================== //
cd "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021"
use "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021\figure1data.dta"

graph bar used2 if urban == 1, over(inc_quintile, sort(inc_quintile_order) label(labsize(large))) b1title("Income Quintile", size(vlarge)) ytitle("Average Diversion Fraction", size(vlarge)) ysc(r(0 0.12)) ylabel(0(0.02)0.12, angle(0) labsize(large)) ytick(0(0.01)0.12) graphregion(color(white)) name(first)

graph bar used2 if urban == 0, over(inc_quintile, sort(inc_quintile_order) label(labsize(large))) b1title("Income Quintile", size(vlarge)) ytitle("Average Diversion Fraction", size(vlarge)) ysc(r(0 0.12)) ylabel(0(0.02)0.12, angle(0) labsize(large)) ytick(0(0.01)0.12) graphregion(color(white)) name(second)



// ============================== Figure 2 ================================== //
clear all

cd "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021"
use "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021\figure2data.dta"

graph bar used2 if car == 1, over(price_group, sort(price_group_order) label(labsize(large))) b1title("Vehicle Price", size(vlarge)) ytitle("Average Diversion Fraction", size(vlarge)) ysc(r(0 0.12)) ylabel(0(0.02)0.12, angle(0) labsize(large)) ytick(0(0.01)0.12) graphregion(color(white)) name(third)

graph bar used2 if car == 0, over(price_group, sort(price_group_order) label(labsize(large))) b1title("Vehicle Price", size(vlarge)) ytitle("Average Diversion Fraction", size(vlarge)) ysc(r(0 0.12)) ylabel(0(0.02)0.12, angle(0) labsize(large)) ytick(0(0.01)0.12) graphregion(color(white)) name(fourth)



// =============================== Figure 3 ================================= //
clear all

cd "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021"
use "L:\Project-MaritzCX\Workspace1\Publication Replication Files\Leard\JAERE 2021\figure3data.dta"

catplot usedage if usedage <= 10, percent vertical b1title("Age of Second Choice Used Vehicle (Years)", size(vlarge)) ytitle("Percent", size(vlarge)) graphregion(color(white)) ylabel(0(5)40, labsize(large)) var1opts(label(labsize(large))) name(fifth)